T=1000;
x=zeros(1,T);
w=randn(1,T);
theta=randn(1)*0.1;
np=length(theta);
P=ones(np);
q=zeros(np,T);
r=ones(1,T);
alpha=1;
for k=np+1:T,
    phi=x(k-p:k-1)';
    x(k)=dot(theta,phi)+w(k);
    r(k)=r(k-1)+phi'*phi;
    alpha=1/log(1+r(k));
    [theta_est,P]=WLS_step(q(:,k-1),P,phi,x(k),alpha);
    q(:,k)=theta_est;
end;
subplot(221);plot(1:T,x);title('x');
subplot(222);plot(1:T,q(1,:));title('\theta_{est}');
subplot(223);plot(1:T,r);title('r');
subplot(224);plot(1:T,w);title('w');
theta
theta_est